#include "AnaDef.h"
#include "style.h"
#include "XGLUtils.h"
#include "AnaUtils.h"

const Double_t LOCALEPSILON = 1e-6;

Int_t ftpc = -999;
Int_t fgrid = -999;
const Int_t fntype = 3;

TList *glist=new TList;

const Int_t fColBl = kBlue;//7; //


Int_t IDx()
{
  const Int_t idx = fgrid/10;
  if(idx!=1 && idx!=2){
    printf("bad idx %d %d\n", fgrid, idx); exit(1);
  }
  return idx;
}

Int_t IDy()
{
  const Int_t idy= fgrid%10;
  if(idy!=0 && idy!=3){
    printf("bad idy %d %d\n", fgrid, idy); exit(1);
  }
  return idy;
}


Double_t getHistX(const Int_t ib, const TH1* hh)
{
  if(IDx()==1){
    return TMath::Sqrt(  hh->GetXaxis()->GetBinLowEdge(ib) * hh->GetXaxis()->GetBinUpEdge(ib) );//TMath::Power(10, (hh->GetXaxis()->GetBinLowEdge(ib)+hh->GetXaxis()->GetBinUpEdge(ib))/2);
    //
  }
  else{
    return hh->GetXaxis()->GetBinCenter(ib);
  }
}

Double_t getHistMinimum(const TH1D *hh)
{
  Double_t ymin = 1e10;
  Int_t bmin = -999;
  for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
    const Double_t yy = hh->GetBinContent(ii);
    if(yy<LOCALEPSILON)
      continue;

    if(yy<ymin){
      ymin=yy;
      bmin = ii;
    }
  }
  if(bmin<0){
    printf("getHistMinimum bmin < 0 \n"); exit(1);
  }
  return ymin;
}

void clearGraph(TGraphAsymmErrors * grs[])
{
  for(Int_t itype=0; itype<fntype; itype++){
    TGraphAsymmErrors* gr = grs[itype];

    //----------------
  const Double_t oldlastx0 = gr->GetX()[gr->GetN()-1];

  Int_t istart = -999;
  for(Int_t ii=0; ii<gr->GetN(); ii++){
    if(gr->GetY()[ii]>LOCALEPSILON){
      istart = ii;
      break;
    }
  }
  if(istart<0){
    printf("clearGraph istart<0\n"); 
    for(Int_t ii = 0; ii < gr->GetN(); ii++){
      printf("raw %d/%d x %.12e %e\n", ii, gr->GetN(), gr->GetX()[ii], gr->GetY()[ii]);
    }
    exit(1);
  }

  const Int_t nn = gr->GetN()-istart;
  for(Int_t ii=0; ii<nn; ii++){
    const Int_t iold = ii+istart;
    const Double_t xx = gr->GetX()[iold];
    const Double_t yy = gr->GetY()[iold];
    const Double_t exl = gr->GetEXlow()[iold];
    const Double_t exh = gr->GetEXhigh()[iold];
    const Double_t eyl = gr->GetEYlow()[iold];
    const Double_t eyh = gr->GetEYhigh()[iold];

    gr->SetPoint(ii, xx, yy);
    gr->SetPointError(ii, exl, exh, eyl, eyh);
  }
  gr->Set(nn);

  const Double_t newlastx0 = gr->GetX()[gr->GetN()-1];
  if(fabs(oldlastx0-newlastx0)>LOCALEPSILON){
    printf("clearGraph old != new %e %e\n", oldlastx0, newlastx0); exit(1);
  }
  //-------------------------
  printf("checkfirstpoint %s %e\n", gr->GetName(), gr->GetX()[0]);
  }
}

Int_t getGrStart(const Double_t x0, const TGraphAsymmErrors * gr)
{
  Int_t start = -999;
  for(Int_t ii = 0; ii < gr->GetN(); ii++){
    if( TMath::Abs(x0-gr->GetX()[ii])< LOCALEPSILON ){
      start = ii;
      break;
    }
  }
  if(start<0){
    printf("getGrStart<0 %d x0 %.12e\n", start, x0); 
    for(Int_t ii = 0; ii < gr->GetN(); ii++){
      printf("raw %d/%d x %.12e\n", ii, gr->GetN(), gr->GetX()[ii]);
    }
    exit(1);
  }
  return start;
}

Int_t getHistStart(const Double_t x0, const TH1* hh)
{
  Int_t start = -999;
  for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
    if( TMath::Abs(x0-getHistX(ii,hh))<LOCALEPSILON  ){
      start = ii;
      break;
    }
  }
  if(start<0){
    printf("getHistStart<0 %d x0 %.12e\n", start, x0);
    for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
      printf("%s %d/%d x %.12e\n",hh->GetName(), ii, hh->GetNbinsX(), getHistX(ii,hh));
    }
    TCanvas * cc = new TCanvas;
    TH1D * htmp = (TH1D*) hh->Clone(Form("testfail%s", hh->GetName()));
    htmp->Draw();
    cc->Print(Form("%s.eps", htmp->GetName()));
    exit(1);
  }
  return start;
}

void multiplyGr(TGraphAsymmErrors * grstat, TGraphAsymmErrors * grsys, TGraphAsymmErrors * tomul, const Bool_t kErr)
{
  if(!kErr){
    printf("bad kerr in multiplyGr %d\n", kErr); exit(1);
  }

  const Double_t x0 = grstat->GetX()[0];

  const Int_t istart = getGrStart(x0, tomul);

  for(Int_t igr=0; igr<grstat->GetN(); igr++){
    const Int_t itom = igr+istart;

    const Double_t xold = grstat->GetX()[igr];
    if(TMath::Abs( xold - tomul->GetX()[itom] )>LOCALEPSILON){
      printf("multiplyGr::grstat x not match! igr %d itom %d stat %e raw %e\n", igr, itom, xold, tomul->GetX()[itom]); exit(1);
    }

    if(TMath::Abs( xold - grsys->GetX()[igr] )>LOCALEPSILON){
      printf("multiplyGr::grsys x not match! igr %d sys %e raw %e\n", igr, grsys->GetX()[igr], tomul->GetX()[itom]); exit(1);
    }
    
    const Double_t yfrac  = grstat->GetY()[igr];
    const Double_t efracl = grstat->GetEYlow()[igr];
    const Double_t efrach = grstat->GetEYhigh()[igr];

    const Double_t ytom  = tomul->GetY()[itom];
    const Double_t etoml = kErr ? tomul->GetEYlow()[itom] : 0;
    const Double_t etomh = kErr ? tomul->GetEYhigh()[itom] : 0;

    //===

    const Double_t ynew = yfrac*ytom;
    const Double_t elnew = ynew*TMath::Sqrt(TMath::Power(efracl/yfrac,2) + TMath::Power(etoml/ytom,2));
    const Double_t ehnew = ynew*TMath::Sqrt(TMath::Power(efrach/yfrac,2) + TMath::Power(etomh/ytom,2));

    //===
    grstat->SetPoint(igr, xold, ynew);
    grstat->SetPointError(igr, grstat->GetEXlow()[igr], grstat->GetEXhigh()[igr], elnew, ehnew);

    //no sys from tomul
    grsys->SetPoint(igr, xold, ynew);
    grsys->SetPointError(igr, grsys->GetEXlow()[igr], grsys->GetEXhigh()[igr], grsys->GetEYlow()[igr]*ytom, grsys->GetEYhigh()[igr]*ytom);
  }
}

void multiplyHist(TGraphAsymmErrors * grstat, TGraphAsymmErrors * grsys, TH1D * tomul, const Bool_t kErr)
{
  if(!kErr){
    printf("bad kerr in multiplyHist %d\n", kErr); exit(1);
  }

  glist->Add(tomul);

  const Double_t x0 = grstat->GetX()[0];

  const Int_t istart = getHistStart(x0, tomul);

  for(Int_t igr=0; igr<grstat->GetN(); igr++){
    const Int_t itom = igr+istart;

    const Double_t xold = grstat->GetX()[igr];
    if(TMath::Abs( xold - getHistX(itom, tomul) )>LOCALEPSILON){
      printf("multiplyGr::grstat x not match! igr %d itom %d stat %e raw %e\n", igr, itom, xold, getHistX(itom, tomul)); exit(1);
    }

    if(TMath::Abs( xold - grsys->GetX()[igr] )>LOCALEPSILON){
      printf("multiplyGr::grsys x not match! igr %d sys %e raw %e\n", igr, grsys->GetX()[igr], xold); exit(1);
    }
    
    const Double_t yfrac  = grstat->GetY()[igr];
    const Double_t efracl = grstat->GetEYlow()[igr];
    const Double_t efrach = grstat->GetEYhigh()[igr];

    const Double_t ytom  = tomul->GetBinContent(itom);
    const Double_t etoml = kErr ? tomul->GetBinError(itom) : 0;
    const Double_t etomh = kErr ? tomul->GetBinError(itom) : 0;

    //===

    const Double_t ynew = yfrac*ytom;
    const Double_t elnew = ynew*TMath::Sqrt(TMath::Power(efracl/yfrac,2) + TMath::Power(etoml/ytom,2));
    const Double_t ehnew = ynew*TMath::Sqrt(TMath::Power(efrach/yfrac,2) + TMath::Power(etomh/ytom,2));

    //===
    grstat->SetPoint(igr, xold, ynew);
    grstat->SetPointError(igr, grstat->GetEXlow()[igr], grstat->GetEXhigh()[igr], elnew, ehnew);

    //no sys from tomul
    grsys->SetPoint(igr, xold, ynew);
    grsys->SetPointError(igr, grsys->GetEXlow()[igr], grsys->GetEXhigh()[igr], grsys->GetEYlow()[igr]*ytom, grsys->GetEYhigh()[igr]*ytom);
  }
}

void addSys(TGraphAsymmErrors * grsys, TH1F * hcorrsys)
{
  glist->Add(hcorrsys);

  TString tmpn(Form("%s%s",grsys->GetName(), hcorrsys->GetName()));
  tmpn.ReplaceAll("sysErrTot","");
  tmpn.ReplaceAll("SysErrTot","");

  TH1F * hrawsys = (TH1F*) hcorrsys->Clone(tmpn+"_rawsys");
  hrawsys->Scale(0);
  glist->Add(hrawsys);

  TH1F * hqsumsys = (TH1F*) hrawsys->Clone(tmpn+"_qsumsys");
  hqsumsys->Scale(0);
  glist->Add(hqsumsys);

  const Double_t x0 = grsys->GetX()[0];

  const Int_t istart = getHistStart(x0, hcorrsys);

  for(Int_t igr=0; igr<grsys->GetN(); igr++){
    const Int_t ics = igr+istart;

    const Double_t xgr = grsys->GetX()[igr];
    if(TMath::Abs( xgr - getHistX(ics, hcorrsys) )>LOCALEPSILON){
      printf("multiplyGr::grsys x not match! igr %d ics %d stat %e raw %e\n", igr, ics, xgr, getHistX(ics, hcorrsys)); exit(1);
    }

    const Double_t yy  = grsys->GetY()[igr];
    const Double_t efracl = grsys->GetEYlow()[igr];
    const Double_t efrach = grsys->GetEYhigh()[igr];

    //yy can be 0 because of low pt CF being cut away and set as 0
    /*
    if(yy<EPSILON){
      printf("addSys strange yy %d %e\n", igr, yy); exit(1);
    }
    */

    if(yy>EPSILON){
      hrawsys->SetBinContent(ics, efracl/yy); hrawsys->SetBinError(ics,0);
    }

    //be careful
    const Double_t ecs = hcorrsys->GetBinContent(ics);

    //===                                                                                                                                                                
    const Double_t elnew = yy<EPSILON? 0 : yy*TMath::Sqrt(TMath::Power(efracl/yy,2) + TMath::Power(ecs,2));
    const Double_t ehnew = yy<EPSILON? 0 : yy*TMath::Sqrt(TMath::Power(efrach/yy,2) + TMath::Power(ecs,2));

    if(yy>EPSILON){
      hqsumsys->SetBinContent(ics, elnew/yy); hqsumsys->SetBinError(ics,0);
    }

    //===                                                                                                                                                                
    grsys->SetPointError(igr, grsys->GetEXlow()[igr], grsys->GetEXhigh()[igr], elnew, ehnew);
  }
}

void DrawHist(const Int_t did)
{
  TCanvas * cc = new TCanvas(Form("tpc%dcid%dgrid%d",ftpc, did, fgrid), "",600,400);
  style::PadSetup(cc);

  gStyle->SetOptTitle(0);
  gPad->SetRightMargin(0.01);
  gPad->SetTopMargin(0.03);

  gPad->SetLogx(0);
  gPad->SetLogy(0);

  Int_t idraw = 0;    
  TLegend * lg = 0x0;

  const Int_t nhist = glist->GetEntries();
  for(Int_t ient = 0; ient<nhist; ient++){
    TH1D *htmp = (TH1D*)glist->At(ient);

    //_______________________________________________________

    const TString hn = htmp->GetName();
    
    if(hn.Contains("ePlusEMinus")){
      continue;
    }
    
    /*
      TCanvas * ctmp = new TCanvas;
      htmp->Draw();
      ctmp->Print(Form("%s.eps", htmp->GetName()));
      delete ctmp;
    */

    //_________________________________________________
    //select correction factor

    TString ytit;
    Double_t hmax =-999;
    Double_t hmin = -999;

    if(did==0 && hn.Contains("BbBCorr")){
      ytit+="bin-by-bin correction factor";
      
      hmax = 3;//21;// 1.75;//21;//3.5;
      hmin = 0;//0.5;//0;//0.5; 
      //gPad->SetLogy(1);

      if(!lg){
        //for global legend
        //lg = new TLegend(0.05, 0.05, 0.95, 0.95);
        lg = new TLegend(0,0,1,1);

        //individual legend
        //lg = new TLegend(0.5, 0.35, 0.95, 0.95);
        //lg = new TLegend(0.35, 0.15, 0.65, 0.55);

        //style::ResetStyle(lg);
        lg-> SetNColumns(3);
        lg->SetFillStyle(-1);
        lg->SetBorderSize(-1);
      }
    }
    else if(did==1 && hn.Contains("SecCorr")){
      ytit+="secondary correction factor";
      
      hmax = 1.05;
      hmin = 0.15;

      if(!lg){
        lg = new TLegend(0.5, 0.20, 0.90, 0.75);
        style::ResetStyle(lg);
      }
    }
    else if(did==2 && hn.Contains("PiMuCorr")){
      ytit+="muon correction factor";
      
      hmax = 1.01;
      hmin = 0.95;

      if(!lg){
        lg = new TLegend(0.5, 0.15, 0.85, 0.55);
        style::ResetStyle(lg);
      }
    }
    else if(did==3 && hn.Contains("ToPiCorr")){
      ytit+="combined correction factor";

      hmax = 2.5;//5;
      hmin = 0;

      if(!lg){
        lg = new TLegend(0.5, 0.5, 0.85, 0.95);
        style::ResetStyle(lg);
      }
    }
    else if(did==4 && (hn.Contains("ysErrTot"))){//"hRelSysErrTot"))){
      ytit+="comb. correction rel. sys. error";

      hmin = 0;
      hmax = 0.3;//0.5;
      if(IDy()==3){
        //hmax = 0.15;
      }

      if(!lg){
        lg = new TLegend(0.5, 0.50, 0.90, 0.95);
        style::ResetStyle(lg);
      }
    }
    else if(did==5 && (hn.Contains("rawsys"))){
      ytit +="raw rel. sys. error";

      hmin = 0;
      hmax = 0.3;

      if(!lg){
        lg = new TLegend(0.5, 0.50, 0.90, 0.95);
        style::ResetStyle(lg);
      }
    }
    else if(did==6 && (hn.Contains("qsumsys"))){
      ytit +="final rel. sys. error";

      hmin = 0;
      hmax = 0.3;

      if(!lg){
        lg = new TLegend(0.5, 0.50, 0.90, 0.95);
        style::ResetStyle(lg);
      }

      //for(Int_t kkk=1; kkk<=htmp->GetNbinsX(); kkk++){
        //printf("testtestkkk %d: %e %e %e\n", kkk, htmp->GetBinCenter(kkk), htmp->GetBinContent(kkk), htmp->GetBinError(kkk));
      //}
    }
    else{
      continue;
    }

    //--
    
    if(hn.Contains("GFL")){
      //ytit+=" (w/ GEANT-FLUKA corr.)";
    }
    
    if(IDy()==0){
      ytit.Prepend("yield ");
    }
    else if(IDy()==3){
      ytit.Prepend("to-pion ratio ");
    }

    //_____________________________________

    style::ResetStyle(htmp);   

    htmp->GetYaxis()->SetTitle(ytit);
    htmp->SetMinimum(hmin);//getHistMinimum(htmp)-htmp->GetBinContent(htmp->GetMaximumBin())*0.1);
    htmp->SetMaximum(hmax);//htmp->GetBinContent(htmp->GetMaximumBin())*1.1);
    
    printf("testtest hn %s did %d\n", hn.Data(), did);
    
    TString xtit;
    if(hn.Contains("Z") || hn.Contains("_z_")){
      xtit="z^{ch}";
      htmp->GetXaxis()->SetRangeUser(0,0.999);
    }
    else if(hn.Contains("Pt") || hn.Contains("_pt_")){
      xtit="p_{T}";
      htmp->GetXaxis()->SetRangeUser(0.15,25);
      htmp->GetXaxis()->SetMoreLogLabels();
      gPad->SetLogx(1);
    }
    else{
      printf("bad hn for xtit %s\n", hn.Data()); exit(1);
    }
    htmp->GetXaxis()->SetTitle(xtit);
    
    htmp->SetMarkerSize(1);

    TString tag;
    Int_t col = kBlack;
    if(hn.Contains("piPlus") || hn.Contains("_pi") || hn.Contains("PiMuCorr")){
      tag+="#pi";
      col = kRed;
    }
    else if(hn.Contains("kPlus") || hn.Contains("_K") || hn.Contains("ktoPi")){
      tag+="K";
      col = fColBl;
    }
    else if(hn.Contains("ppbar") || hn.Contains("_p") || hn.Contains("ptoPi")){
      tag+="p";
      col = kBlack;
    }
    else{
      printf("wrong hn for type %s\n", hn.Data()); exit(1);
    }

    tag+=" ";
    Int_t msty = 24;
    if(hn.Contains("05_10") || hn.Contains("5_10")){
      tag+="p_{T}^{ch jet} 5-10 GeV/c";
      msty = 24;
    }
    else if(hn.Contains("10_15")){
      tag+="p_{T}^{ch jet} 10-15 GeV/c";
      msty = 26;
    }
    else if(hn.Contains("15_20")){
      tag+="p_{T}^{ch jet} 15-20 GeV/c";
      msty = 32;
    }
    else {
      printf("wrong hn for jpt %s\n", hn.Data()); exit(1);
    }

    htmp->SetMarkerStyle(msty);
    htmp->SetMarkerColor(col);
    htmp->SetLineColor(col);

    lg->AddEntry(htmp, tag, "lp");

    if(did>=4){
      htmp->Draw(idraw?"same lp":"lp");
    }
    else{
      htmp->Draw(idraw?"same hist e":"hist e");
    }

    XGLUtils::PrintHist("drawhist", htmp);

    idraw++;
  }

  if(lg&&cc){
    //lg->Draw();
    cc->Print(Form("%s.eps", cc->GetName()));

    if(did==0){
      delete cc;
      cc=new TCanvas("clg","", 650,150);
      style::PadSetup(cc);
      lg->Draw();
      cc->Print(Form("Legend%s.eps", cc->GetName()));
    }
  }

  delete cc;
  delete lg;
}

void ReadRaw(TGraphAsymmErrors * &grraw, const TString jpt)
{
  if(IDy()==0){
    TFile * frawyield = new TFile(Form("rawyield_TPC%d_JPT%s.root", ftpc, jpt.Data()));
    if(!frawyield->IsOpen()){
      printf("no raw yield!! %s\n",  frawyield->GetName()); exit(1);
    }
    frawyield->ls();
    
    grraw = (TGraphAsymmErrors*) frawyield->Get(Form("grid%d", IDx()));
    if(!grraw){
      printf("no grraw!\n"); exit(1);
    }
  }
}

void ReadFrac(TGraphAsymmErrors* grstat[], TGraphAsymmErrors* grsys[], const TString jpt)
{
  TFile * frawfrac = new TFile(Form("rawfraction_TPC%d_JPT%s_GrID%d.root", ftpc, jpt.Data(), fgrid));
  if(!frawfrac->IsOpen()){
    printf("no raw frac!! %s\n",  frawfrac->GetName()); exit(1);
  }
  frawfrac->ls();

  for(Int_t itype=0; itype< fntype; itype++){
    grstat[itype] = (TGraphAsymmErrors*) frawfrac->Get(Form("grstat%d",itype));
    if(!grstat[itype]){
      printf("no grstat %d\n", itype); exit(1);
    }

    grsys[itype] = (TGraphAsymmErrors*) frawfrac->Get(Form("grsys%d",itype));
    if(!grsys[itype]){
      printf("no grsys %d\n", itype); exit(1);
    }

    XGLUtils::PrintGr(Form("readfrac_%s", jpt.Data()),grsys[itype]);
  }
}

//  CheckToPion(hBbB, hSec, hPiMu, hToPi,2);
void CheckToPion(TH1D *hBbB[], TH1D *hSec[], TH1D * hPiMu[], TH1D* hToPi[], const Int_t itype)
{
  Int_t nfatal = 0;
  for(Int_t ibin = 1; ibin<= hToPi[itype]->GetNbinsX(); ibin++){
    const Double_t cfpar  = hBbB[itype]->GetBinContent(ibin) * hSec[itype]->GetBinContent(ibin) *( hPiMu[itype]? hPiMu[itype]->GetBinContent(ibin) : 1);
    const Double_t cfpion = hBbB[1]->GetBinContent(ibin) * hSec[1]->GetBinContent(ibin) *( hPiMu[1]? hPiMu[1]->GetBinContent(ibin) : 1);

    if(fabs( cfpar-cfpion*hToPi[itype]->GetBinContent(ibin)  )>LOCALEPSILON){
      printf("checktopion fail itype %d ibin %d xx %e %e cfpar %e cfpion %e topi %e cfpar/cfpion %e", itype, ibin, hToPi[itype]->GetBinLowEdge(ibin), hToPi[itype]->GetXaxis()->GetBinUpEdge(ibin), cfpar, cfpion, hToPi[itype]->GetBinContent(ibin), cfpar/cfpion); 
      if(hToPi[itype]->GetBinContent(ibin)>LOCALEPSILON){
        printf(" fatal par bbb %e sec %e pimu %e ; pion bbb %e sec %e pimu %e\n", hBbB[itype]->GetBinContent(ibin) , hSec[itype]->GetBinContent(ibin), hPiMu[itype]? hPiMu[itype]->GetBinContent(ibin) : 1, hBbB[1]->GetBinContent(ibin), hSec[1]->GetBinContent(ibin) , hPiMu[1]? hPiMu[1]->GetBinContent(ibin) : 1);
        nfatal++;
      }
      else{
        printf(" still good\n");
      }
    }
  }
  if(nfatal){
    exit(1);
  }
}

void ReadCorrFrac(TH1D *hBbB[], TH1D *hSec[], TH1D * hPiMu[], TH1D* hToPi[], const TString xvar, TString jpt, const Bool_t kSc)
{
  TFile * fcf = new TFile("cf.root");
  if(!fcf->IsOpen()){
    printf("no correction!!\n"); exit(1);
  }
  fcf->ls();

  TDirectory *fdir= fcf->GetDirectory(Form("tpcCut%d",ftpc));
  if(!fdir){
    printf("no fdir!\n"); exit(1);
  }

  const TString pname[]={"ppbarGFL","piPlusPiMinus","kPlusKMinusGFL","ePlusEMinus"};

  jpt.ReplaceAll("5-10","05-10");
  jpt.ReplaceAll("-","_");

  TString stmp;
  const TString sctag=kSc?"":"nonSc_";

  for(Int_t itype=0; itype<fntype; itype++){
    //no Sc tag
    stmp=Form("hBbBCorr%s_%s_%s", xvar.Data(), jpt.Data(), pname[itype].Data());
    hBbB[itype]=(TH1D*) fdir->Get(stmp);
    if(!hBbB[itype]){
      printf("hBbB null!! %d %s\n", itype, stmp.Data());
      fdir->ls();
      exit(1);
    }

    stmp=Form("hSecCorr%s_%s%s_%s", xvar.Data(), sctag.Data(), jpt.Data(), pname[itype].Data());
    hSec[itype]=(TH1D*) fdir->Get(stmp);
    if(!hSec[itype]){
      printf("hSec null!! %d %s\n", itype, stmp.Data());
      fdir->ls();
      exit(1);
    }

    if(itype== AnaUtils::kPION){
      stmp=Form("hPiMuCorr%s_%s%s", xvar.Data(), sctag.Data(), jpt.Data());
      hPiMu[itype]=(TH1D*) fdir->Get(stmp);
      if(!hPiMu[itype]){
        printf("PiMu null %s\n", stmp.Data());
        exit(1);
      }
    }
    else if(itype== AnaUtils::kPROTON || itype== AnaUtils::kKAON){
      stmp=Form("hToPiCorr%s_%s%s_%s", xvar.Data(), sctag.Data(), jpt.Data(), pname[itype].Data()); 
      hToPi[itype]=(TH1D*) fdir->Get(stmp);
      if(!hToPi[itype]){
        printf("no hToPi itype %d %s\n", itype, stmp.Data());
        exit(1);
      }
    }
  }

  //just change to normal scale in case of pt
  if(IDx()==1){
    for(Int_t itype = 0; itype<fntype; itype++){
      TH1D * hhs[]={hBbB[itype], hSec[itype], hPiMu[itype], hToPi[itype]};
      for(Int_t iss=0; iss<(Int_t)(sizeof(hhs)/sizeof(TH1D*)); iss++){
        if(hhs[iss]){
          style::ToNaturalScale(hhs[iss]);
        }
      }
    }
  }

  //void ReadCorrFrac(TH1D *hBbB[], TH1D *hSec[], TH1D * hPiMu[], TH1D* hToPi[], const TString xvar, TString jpt)
  CheckToPion(hBbB, hSec, hPiMu, hToPi, 0);
  CheckToPion(hBbB, hSec, hPiMu, hToPi, 2);
}

void ReadCorrSys(TH1F * hTotCorrSys[], TString xvar, TString jpt)
{
  const TString  pname[]={"p","pi","K"};

  //xvar.ToLower();

  jpt.ReplaceAll("5-10","05-10");
  jpt.ReplaceAll("-","_");

  for(Int_t itype=0; itype<fntype; itype++){
    TString filename;
    TString hname;
    if(IDy()==0){
      filename="outSys_total_species.root";
      hname=Form("hSysErrTot%s_%s_xg_%s", xvar.Data(), jpt.Data(),  pname[itype].Data());
    }
    else if(IDy()==3){
      if(itype==AnaUtils::kPROTON){
        filename="outSys_total_ptoPi.root";
        hname=Form("sysErrTot%s_ptoPi_%s", xvar.Data(), jpt.Data());
      }
      else if(itype==AnaUtils::kKAON){
        filename="outSys_total_ktoPi.root";
        hname=Form("sysErrTot%s_ktoPi_%s", xvar.Data(), jpt.Data());
      }
    }
    if(filename.Length()){
      TFile * fcorrsys = new TFile(filename);
      if(!fcorrsys->IsOpen()){
        printf("no sys correction %s\n", fcorrsys->GetName()); exit(1);
      }
      fcorrsys->ls();
  
      hTotCorrSys[itype]=(TH1F*) fcorrsys->Get(hname);
      if(!hTotCorrSys[itype]){
        printf("no hTotCorrSys itype %d %s\n", itype, hname.Data()); exit(1);
      }

      if(IDx()==1){
        style::ToNaturalScale(hTotCorrSys[itype]);
      }
    }
  }

  if(IDy()==3){
    hTotCorrSys[1] = (TH1F*)hTotCorrSys[0]->Clone("dummy1");
    hTotCorrSys[1]->Scale(0);
  }
  
  if(fntype==4){
    hTotCorrSys[3] = (TH1F*)hTotCorrSys[0]->Clone("dummy3");
    hTotCorrSys[3]->Scale(0);
  }
}

//_______________________________________________________________________________________________________________________________
//_______________________________________________________________________________________________________________________________

void cutoffLowEff(TH1D * hBbB[], TH1D *hSec[], TH1D *hPiMu[], TH1D *hToPi[], TH1F *hTotCorrSys[])
{
  /*
  const TString hn = hbbb->GetName();
  Double_t cutoff = 0.2;
  if(hn.Contains("piPlus") || hn.Contains("ePlus")){
    cutoff = 0.2;
  }
  else if(hn.Contains("kPlus")){
    cutoff = 0.25;
  }
  else if(hn.Contains("ppbar")){
    cutoff = 0.4;
  }
  else if(!hn.Contains("PiMu")){
    printf("wrong hn for cutoff %s\n", hn.Data()); exit(1);
  }

  for(Int_t ii=1; ii<=hbbb->GetNbinsX(); ii++){
    if(hbbb->GetBinLowEdge(ii)<cutoff){
      hbbb->SetBinContent(ii,0);hbbb->SetBinError(ii,0);
    }
  }
  */

  const Double_t cutoffCF = 2;//1.7; //=1/0.6
  const Double_t cutoffCSys = 0.25;
  printf("testtest cutoffCF %e cutoffCSys %e\n", cutoffCF, cutoffCSys);

  for(Int_t itype = 0; itype < fntype; itype++){

    const Int_t nbin = hBbB[itype]->GetNbinsX();  
    Bool_t cutbins[nbin+1];
    
    for(Int_t ibin=1; ibin<=nbin; ibin++){
      cutbins[ibin] = kFALSE;

      if(hBbB[itype]->GetBinContent(ibin) > cutoffCF || hTotCorrSys[itype]->GetBinContent(ibin) > cutoffCSys){
        cutbins[ibin] = kTRUE;
        printf("testtestcutoff  itype %d ibin %d xx %e %e CF %s %e cut %e CS %s %e cut %e\n", itype, ibin, hBbB[itype]->GetXaxis()->GetBinLowEdge(ibin), hBbB[itype]->GetXaxis()->GetBinUpEdge(ibin), hBbB[itype]->GetName(), hBbB[itype]->GetBinContent(ibin), cutoffCF, hTotCorrSys[itype]->GetName(), hTotCorrSys[itype]->GetBinContent(ibin), cutoffCSys);
      }
    }

    /*// bad for hist without CF>2
    //wipe out the startng bins  
    for(Int_t ibin=1; ibin<=nbin; ibin++){
      if(cutbins[ibin])
        break;

      cutbins[ibin] = kTRUE;
    }
    */

    for(Int_t ibin=1; ibin<=nbin; ibin++){
      if(!cutbins[ibin])
        continue;

      const Double_t x0ref = getHistX(ibin, hBbB[itype]);

      TH1 * hhs[]={hBbB[itype],  hSec[itype], hPiMu[itype], hToPi[itype], hTotCorrSys[itype]};
      for(Int_t ih=0; ih<(Int_t)(sizeof(hhs)/sizeof(TH1*)); ih++){
        if(hhs[ih]){

          const Double_t tmpx = getHistX(ibin, hhs[ih]);
          if(fabs(tmpx-x0ref)>LOCALEPSILON){
            printf("bad tmpx %s itype %d ibin %d tmpx %e xref %e\n", hhs[ih]->GetName(), itype, ibin, tmpx, x0ref); exit(1);
          }
          
          printf("testtest %s itype %d ibin %d x0ref %e cut away\n", hhs[ih]->GetName(), itype, ibin, x0ref);
          hhs[ih]->SetBinContent(ibin,0); 
          hhs[ih]->SetBinError(ibin,0); 
        }
      }

    }
  }
}

void combRaw(TGraphAsymmErrors * grraw, TGraphAsymmErrors * grstat, TGraphAsymmErrors * grsys, TH1D * hbbb, TH1D *hsec, TH1D * hpimu)
{
  XGLUtils::PrintGr("posposstat0", grstat);
  XGLUtils::PrintGr("pospossys0", grsys);

  multiplyGr(grstat, grsys, grraw, 1);

  XGLUtils::PrintGr("posposstat1", grstat);
  XGLUtils::PrintGr("pospossys1", grsys);

  multiplyHist(grstat, grsys, hbbb, 1);

  XGLUtils::PrintGr("posposstat2", grstat);
  XGLUtils::PrintGr("pospossys2", grsys);

  multiplyHist(grstat, grsys, hsec, 1);

  XGLUtils::PrintGr("posposstat3", grstat);
  XGLUtils::PrintGr("pospossys3", grsys);

  const TString gn=grstat->GetName();
  if(gn.Contains("1")){//pion
    if(!hpimu){
      printf("muon correction null for pion !!\n"); exit(1);
    }    
    multiplyHist(grstat, grsys, hpimu, 1);
  }
  else if(hpimu){
    printf("strange testtest hpimu %s %s\n", gn.Data(), hpimu->GetName()); exit(1);
  }

  XGLUtils::PrintGr("posposstat4", grstat);
  XGLUtils::PrintGr("pospossys4", grsys);

}

//combToPi(grstat[itype], grsys[itype], hToPi[itype]);
void combToPi(TGraphAsymmErrors * grstat, TGraphAsymmErrors * grsys, TH1D * htopi)
{
  multiplyHist(grstat, grsys, htopi, 1);
}

/*
void CookToPiCS(TH1F *hfracCS[], TH1F *htopiCS[])
{
  for(Int_t itype=0; itype<fntype; itype++){
    if(itype!=AnaUtils::kPROTON && itype!= AnaUtils::kKAON){
      continue;
    }

    htopiCS[itype] = (TH1F*) hfracCS[itype]->Clone(Form("%stopiCS", hfracCS[itype]->GetName()));
    for(Int_t ibin = 1; ibin<=htopiCS[itype]->GetNbinsX(); ibin++){
      const Double_t cspion = hfracCS[AnaUtils::kPION]->GetBinContent(ibin);
      const Double_t cspar  = hfracCS[itype]->GetBinContent(ibin);

      //(1+a)/(1+b)-1 = (a-b)/(1+b), (1-a)/(1-b)-1 = (b-a)/(1-b) --> not symmetric
      //b->0 then |a-b|
      const Double_t vari = TMath::Abs(cspar-cspion);
      htopiCS[itype]->SetBinContent(ibin, vari);
      htopiCS[itype]->SetBinError(ibin, 0);
    }
  }
}
*/

int main(int argc, char* argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=4){
    printf("argc!=3\n"); exit(1);
  }

  ftpc = atoi(argv[1]);

  fgrid = atoi(argv[2]);
  TString xvar;
  if(IDx()==1){
    xvar="Pt";
  }
  else if(IDx()==2){
    xvar="Z";
  }

  const Bool_t kSc = atoi(argv[3]);

  style::SetGlobalStyle(0);

  //_______________________
    
  const TString sjj[]={"5-10","10-15","15-20"};

  for(Int_t ijj = 0; ijj<(Int_t) (sizeof(sjj)/sizeof(TString)); ijj++){
    const TString jpt = sjj[ijj];

    //____________________________________________________________________________________________________________
    //____________________________________________________________________________________________________________
    TGraphAsymmErrors * grstat[]={0x0,0x0,0x0,0x0};
    TGraphAsymmErrors * grsys[]={0x0,0x0,0x0,0x0};
    
    ReadFrac(grstat, grsys, jpt);
    
    //____________________________________________________________________________________________________________
    //____________________________________________________________________________________________________________
    
    TGraphAsymmErrors * grraw = 0x0;
    ReadRaw(grraw, jpt);
    
    //____________________________________________________________________________________________________________
    //____________________________________________________________________________________________________________
    TH1D * hSec[]={0x0,0x0,0x0,0x0};
    TH1D * hBbB[]={0x0,0x0,0x0,0x0};
    TH1D * hPiMu[]={0x0,0x0,0x0,0x0};
    TH1D * hToPi[]={0x0,0x0,0x0,0x0};
    
    ReadCorrFrac(hBbB, hSec, hPiMu, hToPi, xvar, jpt, kSc);
    
    //____________________________________________________________________________________________________________
    //____________________________________________________________________________________________________________
    TH1F * hTotCorrSys[]={0x0,0x0,0x0,0x0}; 

    ReadCorrSys(hTotCorrSys, xvar, jpt);

    //CookToPiCS(hfracCS, htopiCS);

    //____________________________________________________________________________________________________________
    //____________________________________________________________________________________________________________

    //cut away low efficiency for both pt and z
    cutoffLowEff(hBbB, hSec, hPiMu, hToPi, hTotCorrSys);

    //____________________________________________________________________________________________________________
    //____________________________________________________________________________________________________________

    TList * lout = new TList;
    
    for(Int_t itype=0; itype<fntype; itype++){
      lout->Add(grstat[itype]);
      lout->Add(grsys[itype]);
      
      if(IDy()==0){
        combRaw(grraw, grstat[itype], grsys[itype], hBbB[itype], hSec[itype], hPiMu[itype]);
      }
      else if(IDy()==3){
        if(itype!=AnaUtils::kPROTON && itype!=AnaUtils::kKAON){ 
          continue;
        }
        
        combToPi(grstat[itype], grsys[itype], hToPi[itype]);
      }
    
      addSys(grsys[itype], hTotCorrSys[itype]);    
    }

    //____________________________________________________________________________________________________________
    //____________________________________________________________________________________________________________

    XGLUtils::PrintGr("posposstat5",grstat[1]);

    clearGraph(grstat);
    clearGraph(grsys);

    XGLUtils::PrintGr("posposstat5",grstat[1]);

    //+6 for corrected results
    const Int_t idnew = IDy()==0? fgrid+6 : fgrid+4;
    
    //Combined_histRecTHnIDFF_JPT5-10_TPC1TOF9_ETAALL.root 
    TFile * fout =  new TFile(Form("Combined_histRecTHnIDFF_kSc%d_JPT%s_TPC%d_GrID%d.root", kSc, jpt.Data(), ftpc, idnew), "recreate");
    lout->Write();
    fout->Save();
    fout->Close();
    
  }

  //================
  //---

  DrawHist(0);
  DrawHist(1);
  DrawHist(2);
  DrawHist(3);
  DrawHist(4);
  DrawHist(5);
  DrawHist(6);
  //---

  printf("\ncombRawdone!!\n");
  return 0;
}
